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ASYMPTOTICS OF FINITE SYSTEM LYAPUNOV EXPONENTS 
FOR SOME RANDOM MATRIX ENSEMBLES 


PETER J. FORRESTER 


Abstract. For products Pn of N random matrices of size d x d, there is 
a natural notion of finite N Lyapunov exponents {jUi}p=i- In the case of 
standard Gaussian random matrices with real, complex or real quaternion 
elements, and extended to the general variance case for /xi, methods known for 
the computation of linx/v - —>-oo (f^i) are used to compute the large N form of the 
variances of the exponents. Analogous calculations are performed in the case 
that the matrices making up Pn are products of sub-blocks of random unitary 
matrices with Haar measure. Furthermore, we make some remarks relating to 
the coincidence of the Lyapunov exponents and the stability exponents relating 
to the eigenvalues of Pn- 


1. Introduction 


Let Ai, i = 1,..., iV be d x d random matrices independently chosen from the 
same ensemble, and form the product 

Pv = A n A n _i ■ ■ ■ Ai. (1.1) 

Under the condition that the second moment of the diagonal entries of A\Ai are fi¬ 
nite, the multiplicative ergodic theorem of Oseledec [21] [23] tells us that the limiting 
positive definite matrix 


V d := lim (P* P*) 1 /^ 


N—too 

is well defined. In particular Vd has eigenvalues 

e M1 > e M2 > • • • > e w 


( 1 . 2 ) 


(1.3) 


with {ni} referred to as the Lyapunov exponents. 

The recent study [3] , building on [J , gave the asymptotic form of the eigenvalue 
distribution of (P^P/v) 1 ^ 2 ^ in the case that each Ai was chosen from the ensemble 
of d x d standard complex Gaussian matrices. For large N, and with the eigenvalues 
of P(-Pv, {a-’i} say, parametrised by writing Xi = e 2Nyi , this was shown to be given 

by 

j 

(2no?)V* e [LA) 


IfSym n 


N\ 


where, with 'I'(a;) denoting the digamma function 


Mi = ^{d-i + 1), 


= _ 4 >' (d _;+!). 


(1.5) 


In the N —> oo limit, it follows from this and the definitions in the first paragraph 
that the Lyapunov exponents are given by ^{d — i + 1), (i = 1,..., d), which is 
a known result [8j. The value of of for i = 1 agrees with another known result 
[5]. One of the aims of the present paper is to show how {of} for products of 


l 



2 


PETER J. FORRESTER 


Gaussian random matrices can be computed by analogous calculations leading to 
{yo} [IS 0ELZ]. We do this in Section [2721 and we point out there that the mean 
and variance of the finite N Lyapunov exponents for Pn consisting as a product of 
Gaussian and inverse Gaussian matrices follows as a corollary. 

The formulation used to perform the computations in Section 12.21 which itself 
is presented in Section EU allows for the computation of the variances associated 
with the distribution of {z/^} for large but finite N to also be computed in a number 
of other cases. One is for random matrices Ai = E 1 ' 2 ^, where E is a fixed d x d 
positive definite matrix and Gi a standard complex Gaussian random matrix. This 
is carried out in Section 1231 In Section l2~4l we generalise Ai = E 1 / 2 Gi to the case 
that Gi is a standard real, or standard real quaternion, Gaussian random matrix, 
but for the variance corresponding to the largest Lyapunov exponent only. The 
Lyapanov exponents and associated variances are computed for products of sub¬ 
blocks of random unitary matrices with real, complex or real quaternion entries 
and chosen with Haar measure in Section [3] We point out that the formulas for the 
Lyapunov exponents and the associated variances are analogous to those for prod¬ 
ucts involving Gaussian random matrices of two different sizes in equal proportion, 
except for a sign. 

The paper [3] also undertook a study of the eigenvalues of the matrix Pjv, {zk} 
say, parameterised by writing Zk = e NXk+lBk . It was found that for large N the 
joint distribution of (A*,, 9k} is given by (11.411 . with each yi replaced by and mul¬ 
tiplied by l/(27r) d . Thus the angular dependence is uniform, while the asymptotic 
distribution of {A,;} is the same as that for {/Lq}. On this latter point, following [12] 
as reported in the reference book [6, pg. 21], we argue that provided Pn can be 
diagonalised, the Lyaponov exponents can be computed from either the eigenvalues 
of P}jPn or the eigenvalues of Pn itself. This point is made in Section l4~2l 


2. Computation of variances for Gaussian ensembles 

2.1. Known results. While for a general ensemble of d x d random matrices the 
Lyapunov spectrum is difficult to compute, there is a class of probability densities 
on the space of matrices — referred to as isotropic — for which the computation 
becomes relatively simple |5). The defining feature of an isotropic matrix distribu¬ 
tion on { Ai } is that for a particular vector x , Aix/\x\ has distribution independent 
of x. Examples are when Ai = E 1 / 2 c[ l3 ’ d \ where E is a fixed dx d positive definite 
matrix and G^f' d ^ is a d x d standard Gaussian matrix with real (/3 = 1), complex 
(/3 = 2) or real quaternion (/3 = 4) entries. In the latter case, these entries are 
themselves 2x2 matrices with complex entries, and the eigenvalues are calculated 
from the corresponding complex matrix which has a doubly degenerate spectrum; 
see e.g. ;7] §1.3.2]. 

The utility of an isotropic matrix distribution on {Ai} can be seen by revising 
the formalism devised in E3E3] for the computation of the Lyapunov spectrum. 
The main formula specifies the partial sum of Lyapunov exponents according to 

Mi-I-hMfc=sup lim —\ogVol k {yi(N),... ,y k (N)} (k = 1,..., d), (2.1) 

iv—ioo Jy 

where yj{N) := PNUj( 0), the sup operation is over all sets of linearly indepen¬ 
dent vectors {i/i(0),... ,yk{ 0)} and Volt refers to the volume of the parallelogram 
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generated by the given set of k vectors. To quantify the latter, with 

b n ■= [yi(N) y 2 {N) ■ ■ ■ yk{N)}, 

so that J5 jv is the d x k matrix with its columns given by the k vectors {yj{N)}, 
we have 


Vo\ k { yi {N ),..., yk(N)} = det {B^Bn) 1 ' 2 = det^P^P/vBo) 172 - (2-2) 

Suppose now that {yi(0),..., 2 /fc(0)} is a set of k linearly independent unit vec¬ 
tors, and let Edxk denote the d x k matrix with l’s in the diagonal positions of the 
k rows, and 0’s elsewhere. The fact that the distribution of G^f'^x is independent 
of x tells us that 


N 

bZp^PnBo = n E^ xk Gf' d) ^Gf’ d) E dxk 

3 = 1 


N 


=ik 




j,k 


(2.3) 


j =i 


where G^ k ' l> denotes the d x k matrix formed by the first k columns of G^’ d \ and 
it tells us furthermore that the sup operation is redundant. Consequently 

Mr + ■■■+/** = ^ E jf lo & det ( G S d) f EG S d) ) ' ■ ( 2 - 4 ) 

j =i 


while the law of large numbers gives that this can be evaluated as the average over 
{g (M } 

Mi + • • • + Mfc = (log det (G ( J d) t EG ( J d) ) V2 ). (2.5) 

Computation of this average in the case £ = implies union Hz] 

= ^ log | + + !)/ 2 )- (2-6) 

Suppose m is generalised so that each Ai is of size (d + Vi) x (d + Vi-i), uq = 0. 
It has recently been noted in |16 that the calculations leading to (12.611 generalise 
and we have 

N 

Mi = ^ lo g| + ^(^(^-* + 1 )/ 2 )), {V{0j/ 2)) := +j)/2). 

i =1 

(2.7) 

Note that (12.71) is a symmetric function of {u^}. This is in fact a general feature 
of the joint distribution of the eigenvalues of P^Pn [12, and so without loss of 
generality we are free to suppose 0 = u 0 < v\ < V 2 < • • •. The existence of the 
limit limjv^oo u/v implies that the limit in (12.71) is well defined. Hence only a finite 
number of the matrices Ai can be rectangular, with all the rest being square. In 
PS it was furthermore argued, using exact calculations based on knowledge of the 
joint distribution of the spectrum of Pn [21OD1331331 IS! , that the asymptotic form 
of the eigenvalues e 2Nyi in this setting is given by the form ED with fii given by 
(12.71) and 

a 2 = ^(^m-i + l)/2)) ( 2 . 8 ) 

(note that for /3 = 2 — the case of complex entries — this reproduces the formula 
for of in (11.51) 1. Our aim in the next subsection is to show how (12.81) can be deduced 
from analogous calculations as those leading to H1T51 . 
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2.2. Variances for the standard Gaussian ensembles. Our starting point is 
the fact that m is valid for finite N, provided that the limiting operation is 
removed on the RHS, and with {/^} such that the eigenvalues of P^Pn are equal 
to {e 2N ^} (see e.g. the discussion in [3} §6.1]). In the Gaussian case this simplifies 
to (12.41) without the limiting operation on the RHS. The argument leading to (12.51) 
then gives that to leading order in N 

((Ml d -+ Pk) 2 ) — ((pi d - d- Pk)) 2 

= K^( losdet «i ,d)t sG ( j d) ) 1/2 )) - .(logdet ( G r f ^r) 1/2 ) 2 )- 

(2.9) 

The LHS can be rewriten 
k 

Ed+ 2 E {{PiPi) - (Pj)(Pi))- 

j= 1 1 <j<l<k 

But for large IV, and a non degenerate Lyapunov spectrum, we expect the covariance 
( pjpi ) — ( Pj)(pi } to decay exponentially fast. Ignoring such terms, for large N we 
can substitute 

(2-10) 
J= i 

for the LHS of (12.91) . It thus remains to compute the RHS. In the case E = Ifc this 
is a straightforward task. 

Lemma 2.1. We have 

( (log det (G ( J d) f G ( J d) ) V2 ) ) = \(l2 *'((0/2)(d - *) + i/3/2) 

+ ( E( l0 S \ + *(G8/2)(d - *0 + i/3/2)))(2.11) 

i=i ^ 

and thus the RHS of li2.9\) is equal to 

k 

±-^'((ft/2)(d-k)+jft/2). (2.12) 

3 = 1 

Proof. We follow the strategy of the proof given in [8] §2.2] of the evaluation of 
the RHS of (12.51) in the case E = Ifc and j3 = 2. 

The average in (12.111) is over d x k standard Gaussian matrices with real (ft = 1), 
complex (ft = 2) or real quaternion (ft = 4) entries. Such matrices have probability 

_ (R/2') r TV 

density function proportional to e K > k ’ fe , where in the case p = 4 the 

trace operation takes only one the independent eigenvalues (recall the eigenvalues 
in this case are doubly degenerate). In particular, the average is a function of 
(V/Vi) , su gg es ti n g that we change variables by introducing the k x k Wishart 

matrix W^’ k ^ = The corresponding Jacobian is proportional to 

(det w {0 ' k) ) {0/2)(d - k+1 - 2/0) (2.13) 
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(see e.g. [7] eq. (3.22)]). Making use too of the simple identity 

j 2 2 

^(det W^' k) Y ^ = (logdetW (/3 ’ fe) ) 
we see that the LHS of ( 12 . 111 ) is equal to 

i-^/(det 

4 dn 2 V /1 


IV(3A)>0 


(2.14) 


fi=0 


This average is over positive definite Hermitian matrices with real ((3 = 1), complex 
((3 = 2) or real quaternion ((3 = 4) entries, and distribution having a density 
function proportional to 


(det w(/3,fc))(/3/ 2)(d-fc+l-2//3) e -(/3/2)T¥iy(' 3 ' fe ) _ 


(2.15) 


We see that (12.141) is a function only of the eigenvalues of W^ ,k \ suggesting 
that we change variables to the eigenvalues and the eigenvectors. The Jacobian 
which has the feature that the eigenvalue and eigenvector terms factorise — is 
given in e.g. [7] Proposition 1.3.4], and we deduce that (12.141) is equal to 


1 d 2 Z, 


(/3) 


4 dn 2 

r A k,c 0 


fi=0 


(2.16) 


where c = ((3/2)(d — k + fj, + 1 — 2//3) + /r, cq = ((3/2)(d — k + 1 — 2/(3) and 




w , = 

k,a ’ 




dx k e -(fl/ 2 )x ia .a 
1=1 


n 



(2.17) 


The multidimensional integral (12.171) is familiar in random matrix theory as a lim¬ 
iting case of the Selberg integral, and as such is evaluated as a product of gamma 
functions (see e.g. 7J Prop. 4.7.3]). This shows that (12.161) is equal to 


1 _d1 A r(3y» T(((3/2)(d-k) + n + j(3l2) 
4 dM 2 II ^ 2 ) T(((3/2)(d - k) + j(3/2) 


fi—O 


Evaluating the derivative in terms of the digamma function gives (|2.11l) . 

If we now subtract the second term in (12.91) , recalling (12.61) , we see that the term 
in the final line of ( 12 . 111 ) cancels, leaving ( 12 . 121 ) . □ 

Equating (12.101) as the evaluation of the LHS of (12.91) for large N with the 
evaluation (12.121) of the RHS we obtain 


k k 

= - 3 + 1 )/ 2 ), 

i=i 1=1 

which is equivalent to (12.81) in the case of square matrices. 

A minor modification of the above reasoning allows for the derivation of (12.81) 
in the case of rectangular matrices. An important feature is that the non-negative 
integer parameters i/» (* = 1,..., N) can only take on a finite number of different 
values. Let these values be 71 ,..., and suppose they are in proportion ot \ 1 ..., a s 




























6 


PETER .7. FORRESTER 


where Yli=i ch = 1. The RHS of the formula (12.91) with £ = I,j then generalises to 
read 

^ (x>( ( k»g (c<r T ‘' ’G!r +7,) ) I/2 ) 2 ) 

-]>><(( lo S det (of" +7,) *Gf"«'>) 1/2 )) 2 ). (2.18) 
2=1 ' 

Changing variables W^ ,k ^ = Q^ d+ry ^ tgi(/3,d+7>) f or j = 1,... ; s we see that the 
exponent in the Wishart ensemble obtained from the corresponding Jacobian is now 
(/3/2 )(d + 7 i — k + 1 — 2//3). The result of Lemma 12.11 shows that (12.111) holds with 
d replaced by d + u, and all quantities averaged over v in the sense of (1577b . 

A similar averaging of the mean and standard deviation holds in the case that 
Pn consists of both Gaussian and inverse Gaussian matrices (see m for a study 
of the singular values in this setting in the complex case with N finite). From 
the definitions, if the Lyapunov exponents for products of d x d Gaussian random 
matrices are > 4 > '‘‘ > 4 > then the Lyapunov exponents for products 

of d x d inverse Gaussian random matrices are — 4 > ~4-i > ■■■ > — 4- 
Equivalently, denoting the Lyapunov exponents for products of inverse Gaussian 
matrices by /.t)~ > > • • • > , we thus have 

d’i = ~4+\-v (2.19) 

while for the associated variances we have 

K") 2 = (2-20) 

Consider now a matrix product P/v with proportion a + Gaussian matrices and 
a_ inverse Gaussian matrices, and denote the Lyapimov exponents and correspond¬ 
ing variances by {a 4 + ’ - ' ) }, {( cr ,- + _ ' ) ) 2 }- Consideration of the derivation of (j2.18|) . 
and making use of (12.191) and (12.201) . we see that 

/4 + ’ _) = a +4 + a -Vi = a +4 - a-(4 . H _i 

(c r - + ’“ ) ) 2 = a+(cr+) 2 +a_(cr I ") 2 = a+(4? + a -(4+l-if- ( 2 -21) 

2.3. Complex standard Gaussian ensembles with general £. It has been 
shown in 0 that for Ai = £ 1 / 2 Gp ,d ' ) with Gf ‘ ,d ■* a dxd standard complex Gaussian 
random matrix and £ a fixed positive definite matrix with eigenvalues of £ _1 equal 
to {y \,..., yd}, the Lyapunov spectrum corresponding to the product (11.11) is given 

by 

Mfc = - TTTt- ~7 -\ det 

2 I il<i<j<d\y.i Vi) 

This was possible because of the availability of an explicit forrmila for the average 

((det (c3 d) wxyy). 


w 


[(log 

\y) 1 ]i=k+i,-,d 




( 2 . 22 ) 


( 2 . 23 ) 
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Thus, from the works 
e.g. [8| eq. (2.21)]) 

(_l) fc ((fc+l)/ 2 ~d) 


HU we can read off that this average is equal to (see 
[/ 0 °° ^/2+fc-i e -wt dt\ ,=i 




nt-T*'- 


det 


W 


n 


,k 

= 1. d 


i = k + 1,. 
3 = 1... 


(2.24) 


ij =1 (■■■ rii <j<i<d(yj yi ) 

The Harish-Chandra/Itzykson-Zuber matrix integral (see e.g. [7] Prop. 11.6.1]), 
giving a determinant formula for a particular integral over the unitary group is an 
essential component of the derivation of (12.241) . The analogous matrix integral over 
the orthogonal or symplectic group does not permit such a structured evaluation, so 
impeding further progress for calculating the Lyapunov spectrum of the Gaussian 
ensembles with general E and real or real quaternion entries. Differentiating (12.241) 
with respect to /z and setting p, = 0 gives a formula for the right-hand side of (12.91) 
in the case /? = 2. This formula is equivalent to (12.221) . 

For purposes of calculating the variances associated with the Lyapunov spectrum, 
following the strategy of Section 2.2 we require the analogous evaluation of 

2 


logdet(G (2 / )t EG (2 / ) ) 1/2 ) ). 


(2.25) 


Moreover, according to (12.51) the second term in (12.91) is equal to (X^=i(m.?)) 2 ; an d 
we know from (12.221) the value of {)}■ Although the computation of (12.251) is 
elementary, simply requiring a double differentiation of (12.241) with respect to /z, 
unfortunately for general k the resulting formula is not very revealing, consisting of 
a number of terms involving either summations or double summations up to k, with 
some of the summands being determinants. We will therefore restrict attention to 
the cases k = 1, when there is much simplification. 

Proposition 2.2. Let {yi, • • • ,yd} be the eigenvalues of E~, assumed distinct. 
The average \2.‘25\) in the case k = 1 is equal to 


Mi 


T'(l) 1 


d 


(log Vjf 


1 /-A 


4 rL=i,;^ (i-yj/yi) 

It follows that for large N 




log y 3 


- vo/yi) 


(2.26) 


d 


log Vj 


4' , (1)+ / 

4iV V fentW 1 -yj/yi) v ^n Uwii-yj/yi) 


n 


(2.27) 


Proof. Our task is to differentiate (12.241) with k = 1 twice with respect to fa and 

to set y, = 0. Since the only dependence on v for k = 1 is in the first row, this 

reduces to computing the second derivative of each term in the first row. Noting 

d 2 r°° , i / 

— ((log^) 2 - 2(log % )vk(l) + T'(l) + (vk(l)) 2 


— / t^e-^ 1 dt 


d^L 2 


fL—0 


and substituting for the first row of (I2.24[) we see after some simple manipulation 
that the sought average is equal to 


- - - ( 

4 rii<.j<;<d (yi -yj)\ 


det 


[(log 2 /j) 2 ]j=i .rf 

IT 1 *-;-"! 


- 24/(1) det 


[log yj]j=i.d 

[y]- 1 


i=2,...,d 


T'(l) + (T(l)) 2 


(2.28) 


4 


































PETER J. FORRESTER 


where in eliminating the determinant which would otherwise appear on the second 
line use has been made of the Vandermonde determinant evaluation 


de%' = n (*- %')• ( 2 - 2 9) 

1 <j<l<d 

Recalling (12.221) in the case k = 1 allows (12.301) to be equivalently written 


A 4 ! + 


*'( 1 ) 

4 


1 

4ni< j <K d {yi - yj) 


det 


[( lo g%) 2 ]i=i. d 

'.Vj . ■’ 


1 

4 TIi <j<i< d (yi - Vj ) 2 



[log Vj\j=l,...,d 



(2.30) 


Expanding the determinants by the first row and making use of (12.291) to evaluate 
the cofactors we obtain (12.261) . 

To deduce (12.271) from (12.261) we subtract p\ in keeping with (12.91) and equate 
with ( 12 . 101 ) in the case k = 1. □ 


2.4. The variance a\ for Gaussian ensembles with general £. Our deriva¬ 
tion of the formula (12.271) for erf in the case of general variance complex Gaussian 
matrices has made use of the determinant formula (12.241) . which as already men¬ 
tioned has no analogue in the real or real quaternion cases. This was similarly the 
case for the derivation given in [5] of the Lyapunov exponents (12.221) . However, with 
regards to the latter with k = 1 (largest Lyapunov exponent) it has recently been 
shown by Kargin m that an alternative calculation is possible, which furthermore 
does generalise to the real (/? = 1) and real quaternion (/3 = 4) cases, giving the 
formula m Th -!-!] 


2 /ri 


v + 1 ° g (f) + i (*“<“■■> 



dx 

x 


(2.31) 


In (12.311) 7 denotes Euler’s constant and \j = 1 for J true and \J = 0 
Here we adapt the method of m to the computation of the RHS of 


k = 1 . 


otherwise. 
(12.91) with 


Proposition 2.3. Let C denote the simple closed contour starting at large R > 
max{i/i}, goes along the upper edge of the real axis to e > 0 , where e < min{?/i}, 
crosses to the lower edge of the real axis and returns to R. Define 


We have 


JP ~ 27 Ti 


AN V 6 


/ 2 \-0/ 2 dz 

V yi) z 

(2.32) 

■Wi 2 ). 

(2.33) 


Proof. For k = 1 the quantity ^EG'^’^) is a positive scalar random 

variable. Let its density function be denoted pp( A). We then have that 


((logdet(G ( J d)t EG ( J d) ) 1/2 ) 2 ) = \j\\og\fpp{X)d\. 


(2.34) 
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It is shown in m that 




0 -/3Az/2 


II ( Z ~Vi) ' 3/2 ^. Cp = ' 3/2 > ( 2 - 35 ) 


2—1 ' 2=1 

where C is as in (12.321) . Substituting (12.351) in (12.341) . interchanging the order of 
integration and noting 

f(-)~ = 10 * i) 2 -K'O) + -> 2 4 ) 


shows that 


where 


r°° tt 2 

J (log X) 2 pp(\) dX. = ( 7 2 + —)I Q - 27/1 + 1 2 , 


(2.36) 


I p 


1 1 


< 2 - 37 » 

By deforming the contour to go around the pole at 2 = 0 we deduce that Iq = 1. 
We read off from [I7J Eq. (13)] that 

I\ = 2/j.i + 7- (2.38) 

Also, simple manipulation shows 


^2 = (log|) J 0 + 2^1og|)(7i-log|) - Ji 
= (log|) + 2 (log^)( 2 /n+ 7 -log |)- J 2 , 


(2.39) 


where J p is specified by (12.321) . We substitute (12.381) and (12.391) in (12.361) then 
subtract (2/ii) 2 , using the fact that [T7] Eq. (14)] 

2 /ii = -7 + log ~ - Ji (2.40) 

to arrive at (12.331) . □ 

In the case (3 = 2 the integrals in (12.331) can be evaluated using the residue 
theorem to reclaim (12.271) . noting too that ’k'(l) = 7 r 2 / 6 . In the case that y\ = 
• • • = yd = 1 and thus E = I d , provided [3d/2 e Z + use of the residue theorem 
shows that 


dp/ 2-1 1 dp/ 2-1 


2 1 


-Ji= - J i= a * where a * = -Y,-’ ( 2 - 41 ) 

6=1 3 6=2 S p=l P 

which when substituted in (12.331) can be verified to give the same value of cr 2 as 
in the case i = 1 . 

It is shown in P23 that 


Ji = - 


Xa:G( 0 ,l) 


n(‘w) 


x \ ~P/ 2 \ dx 


(2.42) 


as is consistent with the equality between (12.311) and (12.401) . Analogous working 
shows 


J 2 =2 


d 


Xxeco.i) - II ( 4 + f) 

i =1 


X \ ~P/ 2 \ log X Tt 2 

- dx + — 

x 3 


(2.43) 
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so it is also possible to express (12.3311 in terms of real integrals. 

3. Products of random truncated unitary matrices 

3.1. Isotropic property. As remarked in the opening paragraph of EH the 
tractability of the computation of the Lyapunov spectrum for the Gaussian ran¬ 
dom matrices G is due to Gaussian matrices belonging to the class of isotropic 
matrices, meaning that any vector x, Gx/\x\ has distribution independent of x. 
This property holds for any class of random matrices that are unitary invariant, 
and thus their distribution in unchanged by multiplication by a unitary matrix U, 
which is furthermore required to have entries of the same type (real, complex or 
real quaternion) as the random matrices. 

Let Z^ ,d+n ' be a random unitary matrix of size (d + n) x (d + n) with real 
(/3 = 1, complex (j3 = 2) or real quaternion (0 = 4) elements, and chosen with Haar 
measure. The latter has the defining feature that it is unchanged by multiplication 
of the matrices by fixed unitary matrices of the appropriate type. Suppose the fixed 
unitary matrices have the block structure 

rr brfX d 

C — J ! 

A n 

where Vdxd is a d x d unitary matrix. Such matrices act only on the top d x d 
sub-matrix of Z^' d+n \ so it must be that square submatrices of Haar distributed 
unitary matrices (also referred to as truncated unitary matrices [25] ) are themselves 
unitary invariant, and are thus examples of isotropic matrices. With the matrices 
A,j in EH) chosen from the ensemble of truncated unitary matrices, our interest 
in this section is to compute the Lyapunov spectrum as well as the associated 
variances. Crucial to our computations is the fact that a d x k (d > k) sub-block 
Z<-I^ d ) distributed with a probability density function proportional to 

(det(I fc - £ ( J d)t Z ( J d) )) C2 , c 2 := (0/2)(n -k + 1 - 2/0). (3.1) 

We remark that the singular values of products of truncated complex unitary 
matrices for finite N are the topic of the recent work m- 

3.2. Computation of the Lyapunov spectrum and variances. According to 
(12.51) we have 

Hi 3 -f Hk = ( logdet (z { ^ d)] Z { ^ d) ^ 1 (3.2) 

where the average is with respect to m ■ The RHS of (13.21) is easy to compute, 
allowing for a simple formula for the individual /i;. 

Proposition 3.1. Let the matrices A, in (LB be chosen as the top dx d sub-block 
of a (d + n) x (d + n ) Haar distributed random unitary matrix with real (/3 = 1), 
complex (f3 = 2) or real quaternion ((3 = 4) entries. For n> d at least we have 

Hi = \ (*((0/2 )(d - * + 1)) - *((0/2)(n + d - i + 1)) . (3.3) 

Proof. As in the Gaussian case, an essential point is that the average E3 is a 
function only of the positive definite matrix W^’ k ^ = z^ k ' d ^ Z^ k ’ d \ Since Z^ k ’ d ' > is 
a sub-matrix of a unitary matrix, \\Afi-.k) j s fujqhgj- constrained to have eigenvalues 
less than 1, which we denote by writing W^’ k ^ < 1. Regarding this as a change 
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of variables, the corresponding Jacobian is proportional to (12.1311 . Hence our task 
becomes that of evaluating 

i-^-((det w (P,k)y\ , (3.4) 

2dp\ /o<w^’ fc )<i fj,— o 

where has distribution with a probability density function proportional to 

(detW (/3 ’ fc) ) Cl (det(I fc - W (3.5) 

with C 2 as in ED and 

Cl = (/3/2)(d-k + l-2//3). (3.6) 

Both (13.41) and (13.51) are functions of the eigenvalues of W^ ,k ^ only. Introducing 
the appropriate Jacobian as can be found in e.g. [D Proposition 1.3.4] gives for 
(13.41) the expression in terms of multiple integrals 


1 A I (/3) 

-L & 1 fc,Cl+/Z,C2 


2 dn T {f> 

fe,Cl,C2 


(3.7) 


fi =0 


where 


At 3) 

L k,ct 1 ,ct 2 


= / dx i 

Jo 


dx k - xi) c 


i=i 


II \ x l~ x j\ P - ( 3 - 8 ) 

l<j<l<k 

This is precisely the Selberg integral (see e.g. [7J Cli. 4]), already seen in a limiting 
form in (12.171) . Inserting its gamma functions evaluation reduces (13.71) to 


1 d -j—r r(ci + /j. + 1 jp/2)T{ci -T ci + C 2 + 2 + [k + j — l)/3/2) 

2 djl H T(ci + 1 + jp/2)T( Cl + /i + c 2 + 2 + {k + j- 1)0/2) M=0 

k -1 fc-1 

= 2 ( E *( Cl + 1 + 7/3/2) - E *(* + c 2 + 2 + (* + 3 - !)/5/2)). (3.9) 
i=o j=o 

Equating this with the LHS of (13.21) and inserting the values of c\ and C 2 from (13.61) 
and ED we arrive at (1331) . □ 

We remark that even though our derivation of (13.31) requires n > d, used for 
example in the convergence of the integrals in (13.71) . the final expression is well 
defined for n > 0, and we expect it to remain valid. Note in particular that for 
n = 0, (13.31) gives that each Lyapunov is equal to 0. This is consistent with the 
matrices in the product ED then themselves being unitary matrices, as the d x d 
sub-matrix of Z^ ,d+n ' ) is for n = 0 the matrix itself. 

Similar working gives the variances associated with the Lyapunov exponents 

(1331) . 

Proposition 3.2. The analogue of the formula I2.8\) for the variances associated 
with the Lyapunov exponents in the case of products of the truncated unitary ma¬ 
trices of Proposition \3.1\ is given by 

a * = W (^(O 8 / 2 )^ -* + !))- *'((/V2)(n + d-i+ 1)). (3.10) 
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Proof. We use the analogue of the formula (12.911 with the LHS replaced by its 
leading large N form (12.1011 . Thus our main task is to compute 

((logdet (z^'z^y'yy (3.n) 

Following the working in the proof of Proposition 13.11 we see that this average is 
equal to 


1 rl 2 1^ 

k,Ci~\-fl,C2 


4 d/i 2 

k,Cl ,C2 


fi =0 


k -1 


1 d? t— r r(ci + /i + 1 + j/3/2)T(ci + ci + C 2 + 2 + (k + j — l)/3/2) 
4 df-i 2 r(ci + 1 + i/3/2)r(ci + /r + C 2 + 2 + (k + j — l)/3/2) 




/c—1 


fc-1 


i(E^'( Cl + 1 +jP/ 2 ) + C 2 + 2 + {k + j - 

3=0 j =0 

/e—1 fc—1 

+ i(E VE, ( ci + 1 +^/ 2 )-E^ ci + C2 + 2 +( A: +j- 


i=o 


j=o 


We recognise from (13.91) that the final line is equal to (XEi Ai) 2 > an d so cance l s 
out of the RHS of (12.91) . Thus 


k k—1 k —1 

E = 4 ( E *'(* + 1 + W 2 ) - E *'(<* + C2 + 2 + (k + j - 1)0/2)). 

j=i j=o i=o 

Recalling the values of C\ and C 2 from m and m we see that this is equivalent 
to the stated formula (13.101) . □ 

The formulas (13.31) and (13.101) have some similarity with the corresponding for¬ 
mulas (ED and (12.81) for products of rectangular Gaussian random matrices, one of 
size d x d (and thus square), the other of size (d + n) x d. In particular the formulas 
agree if the d x d matrices were to occur in “proportion” 1, and the (d + n) x d 
matrices in “proportion” —1. 


4. Discussion 

4.1. Comparison with numerical experiments. In our earlier paper on Lya¬ 
punov exponents [5] a discussion and references were given on the difficulty of 
accurately computing the Lyapunov exponents. A more robust, but less accurate 
approach is to use Monte Carlo simulation based on m- Numerically stable ways 
to carry out such computations have been reported in [6]. Thus, in relation to /ij, 
and with d = 2, xo = [l 0] , we iteratively define Xi = (i = 1,2,...), 

where y t = Xi/\xi\. According to (12.11) . a Monte Carlo estimate of the largest 
Lyapunov exponent is then given by 
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and similarly an estimate of the corresponding variance is given by 

= ^(^Z]( los l^'l) 2 _ ^i)‘ ( 4 - 2 ) 

3=1 

To use HO to compute p 1 +/ 12 , again taking d = 2 for simplicity, let = [l 0] T 
and Xq 2 - 1 = [0 l] T , then define x- p ^ = Aiy[ p \ (p = 1,2, i = 1,2,...) where 
{ifi iVi^} is obtained from by the Gram-Schmidt orthonormalization 

procedure. It follows from ED that a Monte Carlo estimate of /ii + /12 is 

1 N 

Mi + 1*2 = — ^logdet(XjXj) 1/2 , (4.3) 

3=1 

where A/ is the 2x2 matrix with columns given by x ^ and x[ 2 \ while an estimate 
of the corresponding variance is 

^ 1+^2 = ^(^5Z( logdet ^) X j) 1/2 ) 2 “ (Mi +M2) 2 )- (4.4) 

3 =1 


As an example, let’s consider the case of (11.11) with A,; = E 1 / 2 G[ 2,2 \ We choose 
the eigenvalues of E _1 to equal 1,1/4 so that 


A: = 


1 0 
0 2 


G 


( 2 , 2 ) 


The formula (12.271) gives Naf = — |(log2) 2 ^ = 0.197699..., whereas a Monte 

Carlo calculation based on (14.21) with N = 10 6 gave Naf ss 0.19686. 

To perform Monte Carlo simulations in relation to products of truncated unitary 
matrices, we must first generate random unitary matrices with Haar measure. One 
way to do this is to apply the Gram-Schmidt orthogonalization producure to a 
standard Gaussian matrix. In the complex case, with d = n = 2 we did this and 
with N = 10 5 found the estimates pi ~ —0.41672 and Mi + M 2 ~ —1.6685 whereas 
according to ProDOsition l3.il /ii = — ^ = —0.4166... and pi +/i 2 = —| = —1.666.... 
For the corresponding variances we found Na\ ~ 0.0894 and iV(cr 2 + cr|) ~ 0.3952 
whereas Proposition 13.21 gives Naf = 13/144 = 0.0902... and IV(<t 2 -t-crf) = 29/72 » 
0.4027.... 


4.2. Asymptotic spectrum of Pjy. As revised in the first paragraph of the In¬ 
troduction, the Lyapunov exponents are defined in terms of the eigenvalues of the 
positive definite matrix PjfPjy, where Pn is the random matrix product (11.11) . In 
the recent works [3 HU, for the particular cases that the A,; in m are standard 
Gaussian matrices with real, complex or real quaternion entries, explicit asymp¬ 
totic analysis of the joint distribution of the spectrum of Pjy has revealed that with 
the eigenvalues written in the form Zk = e NXk+l6k , the {A&} — referred to as the 
stability exponents — are identical to the Lyapunov exponents. Moreover in [3] 
a proof of the equality between the Lyapunov and stability exponents in the case 
of 2 x 2 random matrices from an isotropic ensemble has been given. It was then 
conjectured that all matrix products formed from an isotropic ensemble have this 
property. 
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Actually the question of the asymptotic equality between the Lyapunov and sta¬ 
bility exponents has been investigated at least as far back as 1987 [12], albeit in 
the setting of the linearization of sets of differential equations about fixed points. 
In this work it was conjectured that the exponents agree whenever the spectrum of 
Pn is non-degenerate. The reference book ]6| cites [12], and furthermore addresses 
directly the case of products of random matrices, writing “If the matrix Pjv can be 
written in diagonal form, then [in our notation] A*, = /i/A. One immediate conse¬ 
quence of this conjecture is the prediction of the Lyapunov spectrum for a product 
of upper triangular random matrices. There, up to ordering {A*,} = {log(| Afcfc ]}} 
where Akk denotes the entry ( kk ) in the matrices A, of (11.11) . so the conjecture 
gives this same set for the Lyapunov exponents. In the case of /ii this agrees with 
a known theorem [22]. 

Following [13] . we emphasize that the conjecture A k = Hk appears to be an 
intrinsic property of the asymptotic limit of random matrix products, and is not 
shared by other asymptotic limits of spectra in random matrix theory. For example, 
in the case of N x N standard real, complex or real quaternion standard Gaussian 
random matrices 

lim = 

iV—>• OO 

where vi{Y) denotes the eigenvalue of largest modulus of Y. That this ratio is at 
least unity follows from the general inequality \v\(Y)\ < v\ ((F^F) 1 / 2 ), which in 
turn is a consequence of the sub-multiplicity of the RHS viewed as a matrix norm. 

Acknowledgements. This work was supported by the Australian Research Coun¬ 
cil for the project DP140102613. I thank Jesper Ipsen for comments on the first 
draft. 
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